Poincare oscilator¶

$$ \begin{equation} \dot{x} = \gamma * (A - r) * x - \omega * y \end{equation} $$ $$ \begin{equation} \dot{y} = \gamma * (A - r) * y + \omega * x \end{equation} $$ $$ \begin{equation} r = \sqrt{x^2 + y^2} \end{equation} $$

In [1]:
import sympy as sp

# Define the variables and parameters
x, y, gamma, A, v = sp.symbols('x y gamma A v', real=True)
r = sp.sqrt(x**2 + y**2)

# Define r in terms of x and y
r = sp.sqrt(x**2 + y**2)

# Define the system of equations
dxdt = gamma * (A - r) * x - 2 * sp.pi * v * y
dydt = gamma * (A - r) * y + 2 * sp.pi * v * x
In [2]:
# Compute the Jacobian matrix
J = sp.Matrix([[dxdt.diff(x), dxdt.diff(y)],
               [dydt.diff(x), dydt.diff(y)]])
J
Out[2]:
$\displaystyle \left[\begin{matrix}- \frac{\gamma x^{2}}{\sqrt{x^{2} + y^{2}}} + \gamma \left(A - \sqrt{x^{2} + y^{2}}\right) & - \frac{\gamma x y}{\sqrt{x^{2} + y^{2}}} - 2 \pi v\\- \frac{\gamma x y}{\sqrt{x^{2} + y^{2}}} + 2 \pi v & - \frac{\gamma y^{2}}{\sqrt{x^{2} + y^{2}}} + \gamma \left(A - \sqrt{x^{2} + y^{2}}\right)\end{matrix}\right]$

Sled in lastne vrednosti¶

$$ J = \begin{bmatrix} -\frac{\gamma x^2}{r} + \gamma (A - r) & -\gamma \frac{x y}{r} - \omega \\ -\gamma \frac{x y}{r} + \omega & -\frac{\gamma y^2}{r} + \gamma (A - r) \end{bmatrix} $$

$$ Tr(J) = 2 \gamma A - 3 \gamma r = \gamma (2 A - 3 r) $$

Fiksne točke¶

$$ \begin{bmatrix} \dot{x}\\ \dot{y} \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix} = \begin{bmatrix} \gamma * (A - r) * x - \omega * y\\ \gamma * (A - r) * y - \omega * x \end{bmatrix} $$

$$ \omega * y = \gamma * (A - r) * x $$

$$ -\omega * x = \gamma * (A - r) * y $$

Dobimo dve rešitvi

  1. $x = y = 0$

  2. $ r = A $

Stabilnost¶

Točka $x_1 = y_1 = 0$ je stabilna, če je $Tr(J) > 0$ in $Det(J) > 0$ $$ Tr(J)|_{x1,y1} = 2*\gamma*A $$ $$ Det(J)|_{x1,y1} = \gamma^2 A^2 + \omega^2 $$

Vidimo lahko, da velja Tr $Tr(J) > 0$ in $Det(J) > 0$. To pomeni, da je točka $x_1 = y_1 = 0$ nestabilna.

Točka $r = A$. $$ Tr(J)|_{x1,y1} = -\gamma*A $$ $$ Det(J)|_{x1,y1} = \omega^2 $$

Vidimo lahko, da velja Tr $Tr(J) > 0$ in $Det(J) > 0$. To pomeni, da je točka $r = A$ nestabilna.

In [3]:
# Substitute r = A (to evaluate on the limit cycle where amplitude is constant)
J_on_limit_cycle = J.subs(x, 0).subs(y,0)
print("\nJacobian on the limit cycle (x = y = 0):")
J_on_limit_cycle
Jacobian on the limit cycle (x = y = 0):
Out[3]:
$\displaystyle \left[\begin{matrix}A \gamma & - 2 \pi v\\2 \pi v & \text{NaN}\end{matrix}\right]$
In [4]:
# Substitute r = A (to evaluate on the limit cycle where amplitude is constant)
J_on_limit_cycle = J.subs(r, A)
print("\nJacobian on the limit cycle (r = A):")
J_on_limit_cycle
Jacobian on the limit cycle (r = A):
Out[4]:
$\displaystyle \left[\begin{matrix}- \frac{\gamma x^{2}}{A} & - 2 \pi v - \frac{\gamma x y}{A}\\2 \pi v - \frac{\gamma x y}{A} & - \frac{\gamma y^{2}}{A}\end{matrix}\right]$
In [5]:
eigenvalues = J_on_limit_cycle.eigenvals()
print("\nEigenvalues of the Jacobian on the limit cycle (r = A):")
eigenvalues
Eigenvalues of the Jacobian on the limit cycle (r = A):
Out[5]:
{-gamma*(x**2 + y**2)/(2*A) - sqrt((-4*pi*A*v + gamma*x**2 + gamma*y**2)*(4*pi*A*v + gamma*x**2 + gamma*y**2))/(2*A): 1,
 -gamma*(x**2 + y**2)/(2*A) + sqrt((-4*pi*A*v + gamma*x**2 + gamma*y**2)*(4*pi*A*v + gamma*x**2 + gamma*y**2))/(2*A): 1}
In [6]:
trace = J_on_limit_cycle.trace()
determinant = J_on_limit_cycle.det()
print("\nTrace:")
trace
Trace:
Out[6]:
$\displaystyle - \frac{\gamma x^{2}}{A} - \frac{\gamma y^{2}}{A}$
In [7]:
#Determinant of J
print("\nDeterminant:")
determinant
Determinant:
Out[7]:
$\displaystyle 4 \pi^{2} v^{2}$
In [ ]:
import numpy as np
import matplotlib.pyplot as plt

# Parameters
A = 1.0  # Example value for A
gamma = 1.0  # Example value for gamma

# Define a grid of x and y values
x_vals = np.linspace(-2, 2, 100)
y_vals = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x_vals, y_vals)

# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
No description has been provided for this image
In [9]:
import numpy as np
import matplotlib.pyplot as plt
class Oscillator:
    def __init__(self, index, gamma, A, v, initial_x=0.0, initial_y=0.0):
        self.index = index
        self.gamma = gamma
        self.A = A
        self.v = v
        self.x = initial_x
        self.y = initial_y

    def compute_radius(self):
        return np.sqrt(self.x**2 + self.y**2)

    def update_dynamics(self, coupling_x, coupling_y, epsilon, dt):
        r = self.compute_radius()

        # Compute derivatives based on the equations
        dxdt = self.gamma * (self.A - r) * self.x - 2 * np.pi * self.v * self.y + epsilon * coupling_x
        dydt = self.gamma * (self.A - r) * self.y + 2 * np.pi * self.v * self.x + epsilon * coupling_y

        # Update the state using Euler's method
        self.x += dxdt * dt
        self.y += dydt * dt

    def get_phase(self):
        return np.arctan2(self.y, self.x)

    def get_state(self):
        return self.x, self.y
In [10]:
# Initialize an oscillator at position (i, j) with specific parameters
gamma = 10.0
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()
No description has been provided for this image
In [28]:
# Initialize an oscillator at position (i, j) with specific parameters
gamma = 1.0
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()
No description has been provided for this image
In [11]:
# Initialize an oscillator at position (i, j) with specific parameters
gamma = 0.5
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()
No description has been provided for this image
In [12]:
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm.auto import tqdm as tqdm

# Parameters
N = 40  # Size of the lattice (NxN)
gamma = 1.0  # Self-regulation rate
A = 1.0  # Desired amplitude of oscillations
omega = 1.0  # Natural frequency of oscillators
epsilon = 0.1  # Coupling strength
dt = 0.01  # Time step for integration
steps = 2000  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function
def simulate_lattice():
    phases = np.zeros((N, N, steps))
    
    for t in tqdm(range(steps),desc='Simulating lattice'):
        # Update each oscillator
        for i in range(N):
            for j in range(N):
                # Calculate coupling from neighbors
                neighbors = [
                    lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                    lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
                ]
                coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
                coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
                
                # Update oscillator dynamics with coupling
                lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
                
                # Record the phase
                phases[i, j, t] = lattice[i, j].get_phase()
    
    return phases

# Run the simulation
phases_over_time = simulate_lattice()

# Set up the figure and axis for the animation
fig, ax = plt.subplots()
im = ax.imshow(np.sin(phases_over_time[:, :, 0]), cmap="twilight", origin="lower")
plt.colorbar(im, ax=ax, label="Phase (sin)")

# Animation function
def update(frame):
    im.set_array(np.sin(phases_over_time[:, :, frame]))
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [im]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# To display the animation inline (if using Jupyter Notebook), use:
# from IPython.display import HTML
# HTML(ani.to_jshtml())

# To save the animation as a .gif or .mp4 file, uncomment one of these:
# ani.save("poincare_oscillators.gif", writer="imagemagick")
# ani.save("poincare_oscillators.mp4", writer="ffmpeg")

plt.show()
Simulating lattice:   0%|          | 0/2000 [00:00<?, ?it/s]
No description has been provided for this image
In [ ]:
from IPython.display import HTML
HTML(ani.to_jshtml())
Animation size has reached 20977680 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.
Out[ ]:
No description has been provided for this image
The Kernel crashed while executing code in the current cell or a previous cell. 

Please review the code in the cell(s) to identify a possible cause of the failure. 

Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. 

View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.
In [36]:
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
N = 20  # Size of the lattice (NxN)
gamma = 1.0  # Self-regulation rate
A = 1.0  # Desired amplitude of oscillations
omega = 1.0  # Natural frequency of oscillators
epsilon = 0.1  # Coupling strength
dt = 0.01  # Time step for integration
steps = 500  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function
def simulate_lattice():
    phases = np.zeros((N, N, steps))
    
    for t in range(steps):
        # Update each oscillator
        for i in range(N):
            for j in range(N):
                # Calculate coupling from neighbors
                neighbors = [
                    lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                    lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
                ]
                coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
                coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
                
                # Update oscillator dynamics with coupling
                lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
                
                # Record the phase
                phases[i, j, t] = lattice[i, j].get_phase()
    
    return phases

# Run the simulation
phases_over_time = simulate_lattice()

# Set up the figure and axis for the animation
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
scatter = ax.scatter([], [], s=30, c='blue', alpha=0.7)  # Scatter points for phases
arrow = ax.annotate('', xy=(0, 0), xytext=(0, 0), arrowprops=dict(facecolor='red', shrink=0.05, width=2))

# Limit the plot to the unit circle
ax.set_ylim(0, 1)

# Animation update function
def update(frame):
    phases = phases_over_time[:, :, frame].ravel()
    x = np.cos(phases)
    y = np.sin(phases)
    
    # Update scatter points in polar coordinates
    scatter.set_offsets(np.c_[phases, np.ones_like(phases)])  # (angle, radius)

    # Calculate the order parameter (mean phase vector)
    order_param = np.mean(np.exp(1j * phases))  # Mean of complex phases
    order_angle = np.angle(order_param)         # Mean phase angle
    order_magnitude = np.abs(order_param)       # Coherence measure (0 to 1)

    # Update the arrow (order parameter vector)
    arrow.set_position((order_angle, order_magnitude))
    arrow.xy = (order_angle, order_magnitude)
    
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [scatter, arrow]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# To display the animation inline (if using Jupyter Notebook), use:
# from IPython.display import HTML
# HTML(ani.to_jshtml())

# To save the animation as a .gif or .mp4 file, uncomment one of these:
# ani.save("poincare_oscillators_polar.gif", writer="imagemagick")
# ani.save("poincare_oscillators_polar.mp4", writer="ffmpeg")

plt.show()
No description has been provided for this image
In [37]:
# To display the animation inline (if using Jupyter Notebook), use:
from IPython.display import HTML
HTML(ani.to_jshtml())
Animation size has reached 21024736 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.
Out[37]:
No description has been provided for this image
In [39]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from ipywidgets import interact, FloatSlider
from IPython.display import HTML

class Oscillator:
    def __init__(self, index, gamma, A, v, initial_x=0.0, initial_y=0.0):
        self.index = index
        self.gamma = gamma
        self.A = A
        self.v = v
        self.x = initial_x
        self.y = initial_y

    def compute_radius(self):
        return np.sqrt(self.x**2 + self.y**2)

    def update_dynamics(self, coupling_x, coupling_y, epsilon, dt):
        r = self.compute_radius()
        dxdt = self.gamma * (self.A - r) * self.x - 2 * np.pi * self.v * self.y + epsilon * coupling_x
        dydt = self.gamma * (self.A - r) * self.y + 2 * np.pi * self.v * self.x + epsilon * coupling_y
        self.x += dxdt * dt
        self.y += dydt * dt

    def get_phase(self):
        return np.arctan2(self.y, self.x)

    def get_state(self):
        return self.x, self.y

# Parameters
N = 20  # Size of the lattice (NxN)
gamma = 1.0  # Initial self-regulation rate
A = 1.0  # Initial desired amplitude of oscillations
omega = 1.0  # Initial natural frequency of oscillators
epsilon = 0.1  # Initial coupling strength
dt = 0.01  # Time step for integration
steps = 500  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function for a single time step
def simulate_step():
    for i in range(N):
        for j in range(N):
            # Calculate coupling from neighbors
            neighbors = [
                lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
            ]
            coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
            coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
            
            # Update oscillator dynamics with coupling
            lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
    
    # Return the current phases for visualization
    return np.array([[osc.get_phase() for osc in row] for row in lattice])

# Set up the figure and axis for the animation
fig, ax = plt.subplots()
im = ax.imshow(np.sin(simulate_step()), cmap="twilight", origin="lower")
plt.colorbar(im, ax=ax, label="Phase (sin)")

# Animation update function
def update(frame):
    phases = simulate_step()  # Run one step of the simulation with current parameters
    im.set_array(np.sin(phases))  # Update the plot with the new phase values
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [im]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# Display the animation inline in Jupyter
HTML(ani.to_jshtml())

# Interactive sliders using ipywidgets for real-time parameter adjustment
@interact(gamma=FloatSlider(value=gamma, min=0.0, max=2.0, step=0.1),
          A=FloatSlider(value=A, min=0.5, max=2.0, step=0.1),
          omega=FloatSlider(value=omega, min=0.5, max=2.0, step=0.1),
          epsilon=FloatSlider(value=epsilon, min=0.0, max=1.0, step=0.05))
def update_params(gamma, A, omega, epsilon):
    # Update the global variables for gamma, A, omega, and epsilon
    for row in lattice:
        for osc in row:
            osc.gamma = gamma
            osc.A = A
            osc.v = omega
    # This function doesn't return anything but will refresh the parameters.
No description has been provided for this image
interactive(children=(FloatSlider(value=1.0, description='gamma', max=2.0), FloatSlider(value=1.0, description…
In [ ]: